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Abstract. The Surface Cauchy-Born (SCB) method is a computational multi-scale 
method for the simulation of surface-dominated crystalline materials. We present an 
error analysis of the SCB method, focused on the role of surface relaxation. 

In a linearized ID model we show that the error committed by the SCB method is 
0(1) in the mesh size; however, we are able to identify an alternative "approximation 
parameter" — the stiffness of the interaction potential — with respect to which the 
error in the mean strain is exponentially small. Our analysis naturally suggests an 
improvement of the SCB model by enforcing atomistic mesh spacing in the normal 
direction at the free boundary. 



1. Introduction 

Miniaturization of materials to the nanometer scale has led to unexpected and often 
enhanced mechanical properties that are not found in corresponding bulk materials 
jH [28] . The size-dependence of the mechanical behavior and properties has been ex- 
perimentally observed to begin around a scale of about 100 nanometers [18]. A fully 
atomistic simulation of a nanostructure of this size would require on the order of 10 8 
atoms, which motivates the need for computationally efficient multiscale methods. 

The underlying cause for the size-dependent mechanical properties is that surface 
atoms have fewer bonding neighbours, or a coordination number reduction, as com- 
pared to atoms that lie within the material bulk. This results in the elastic properties 
of surfaces being different from those of an idealized bulk material [18] , which becomes 
important with decreasing structural size and increasing surface area to volume ratio [I] . 
Additionally, nanoscale surface stresses [3], which also arise from the coordination num- 
ber reduction of surface atoms [30] , cause deformation of not only the surfaces, but also 
the underlying bulk [H], and can result in unique physical properties such as phase 
transformations [5], or shape memory and pseudoelasticity effects in FCC nanowires 
that are not observed in the corresponding bulk material [35] . 

To study surface-dominated nanostructures, Park et al. recently developed the sur- 
face Cauchy-Born (SCB) model [Ml EH G2]- The idea is to seek an energy functional 



Date: December 6, 2011. 

2000 Mathematics Subject Classification. 70C20, 70-08, 65N12, 65N30. 

Key words and phrases, surface-dominated materials, surface Cauchy-Born rule, coarse-graining. 

KJ and CM were supported by undergraduate vacation bursaries at the Oxford Centre for Nonlinear 
PDE. CO was supported by the EPSRC Grant EP/H003096 "Analysis of Atomistic-to-Continuum 
Coupling Methods". HP was supported by NSF grants CMMI-0750395 and CMMI-1036460. 



2 K. JAYAWARDANA, C. MORDACQ, C. ORTNER, AND H. S. PARK 




5 10 15 20 25 30 5 10 15 20 25 30 

I t 



Figure 1. Displacements and displacement gradients of an atomistic so- 
lution and a surface Cauchy-Born solution, relative to the bulk Cauchy- 
Born solution, for a ID model problem. We observe unexpectedly high 
accuracy at the finite element nodes despite a large error in the displace- 
ment gradient. 



of the form 

E sch (y)= [ W(dy)dx + [ ^(dy,u)ds, 
Jn Jan 

where Q C M 3 is an elastic body, y : Q — > M 3 a deformation field, W the bulk stored 
energy function, and 7 a surface stored energy function. The potentials W, 7 are chosen 
such that VF(F) denotes the energy per unit volume in an infinite crystal under the 
deformation y(x) = Fx, while 7(F, u) is the surface energy per unit area of a half-space 
with surface normal v, under the deformation y(x) = Fx. Thus, W and 7 are derived 
from the underlying atomistic model. For W this is a well-understood idea [TJ [8] ; the 
novel approach in the SCB method is to apply the same principle to the surface energy 
potential. 

In contrast to the SCB method, most computational models (see, e.g., [TUJ, ) are 
based upon a finite element discretization of the governing surface elasticity equations 
of Gurtin and Murdoch [UJ, where the constitutive relation for the surface is linearly 
elastic or uses standard hyperelastic strain energy functions [13J. 

The SCB model was successfully applied to various nanomechanical boundary value 
problems, including thermomechanical coupling [33], resonant frequencies, and elucidat- 
ing the importance of nonlinear, finite deformation kinematics on the resonant frequen- 
cies of both FCC metal [23] and silicon nanowires [TH] IT?] , bending of FCC metal [34J 
nanowires, and electromechanical coupling in surface-dominated nanostructures [T5] . 

The purpose of the present work is to initiate a mathematical analysis of the accuracy 
of the SCB method. We focus on the simplest setting where the only effect is a surface 
relaxation in normal direction. While the SCB model does include surface physics that 
are neglected in the standard Cauchy-Born (CB) model, due to employing a coarse 
finite element discretisation it does not resolve the resulting boundary layer; see the 
numerical results in [9] as well as see Figure [l] for a ID toy model demonstrating this. 
It is therefore a priori unclear to what extent the SCB improves upon the CB model. 
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Figure [I] suggests that, while the error in the displacement and displacement gradient 
is indeed of order 0(1) in the boundary layer, the displacement error in finite element 
nodes is visually negligable, which would imply that the SCB model approximates 
the mean strain (and possibly other averaged quantities) to a much higher degree of 
accuracy. This was indeed observed in extensive numerical tests presented in [2T1 1221 19] . 

There is no traditional discretisation or approximation parameter available with re- 
spect to which we might try to explain this effect. Instead, our analysis measures the 
SCB error in terms of the stiffness of the interaction potential. This enables us to 
identify a suitable asymptotic limit for our analysis on a linearized model problem. We 
confirm the analytical predictions with numerical experiments on the fully nonlinear 
problem in ID and and a periodic semi-infinite 2D domain. 

To the best of our knowledge, our work presents the first approximation error results 
for the SCB method. Although our analysis is elementary, it makes two important novel 
contributions: 1. We show that the "correct" approximation parameter is the stiffness 
of the interaction potential (however, Theil [32J uses similar ideas for an analysis of 
surface relaxation); and 2. We show that the mean strain (which is an important 
quantity of interest) has a much lower error than the strain field. 3. Our results 
show how to substantially improve the accuracy of the SCB method at little additional 
computational cost. Finally, we hope that this work will stimulate further research on 
computationally efficient multiscale methods for surface-dominated nanostructures. 

The issues we address here are closely related to the classical problem of numerical 
methods for resolving boundary layers [25]. The main difference in our case is the 
discrete setting which does not give us the opportunity to let the mesh-size tend to 
zero. For a mathematical analysis of thin atomistic structures, surface energies and 
surface relaxation we refer to [271 121 |32j [26] and references therein. Our work also draws 
inspiration from [Ul [7J where a similar linearised model problem is used to analyze the 
accuracy of atomistic-to-continuum coupling methods. 



2. Analysis of a ID Model Problem 

2.1. Atomistic model. We consider a semi-infinite chain of atoms with reference po- 
sitions £ G N, and deformed positions i/f, f 6 N. We assume that the chain interacts 
through second-neighbour Morse pair interaction. Hence, a deformed configuration y 
has energy 

oo 

where is a shifted Morse potential with stiffness parameter a > and potential 
minimum r$ > 0, 

(j){r) = exp(— 2a(r — r )) — 2exp(— a(r — r )) — (f>o, 

where O is chosen to that W(l) = 0, where W(r) := 4>{r) + 0(2r) = 0, r is defined 
such that W'{1) = 0, 

1 / 1 + 2e~ a \ 
r » = 1 + a'° g (lT2^)- (2) 
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and a > 1 + remains a free parameter. This restriction on a ensures that <f>"(2) < 0, 
which will be convenient in the analysis. The shift of the potential by 0o ensures that 
E a is well-defined. 

The potential W is called the Cauchy-Born stored energy density. We have chosen 
the parameters in the Morse potential so that 1 is the minimizer of W, that is, we are 
working in non-dimensional atomic units. 

Since E a is translation invariant, it is convenient to fix yo = 0. In that case, yi is 
completely determined by the forward differences y' e := ye+\ — ye. Hence we change 
coordinates from the deformation ye to the displacement gradient ut :— y[ — 1, and 
rewrite E & as 

oo 

E a (u) := [<K X + u i) + <t>( 2 + u e + ^+i) • 

1=0 

The proof of the next result, which establishes that E & is well-defined, is given in the 
appendix. 

Proposition 1. E a is well-defined and twice Frechet differentiate in £ 1 (N) with first 
and second variations given by 

oo 

(5E a (u),v) = [<P'(l + u e )v e + <p'{2 + ue + ue +1 ){ve + v e+1 )j, 

oo 

(5 2 E a (u)v, w) = ^ + u t) v e w i + 0"( 2 + u £ + u e+ i)(v e + v i+1 ){w t + w i+ i) ■ 
e=o 



2.2. The Cauchy— Born and surface Cauchy— Born models. The Cauchy Born 
approximation is designed to model elastic bulk behaviour in crystals. The stored 
energy density is chosen so that the Cauchy-Born energy is exact under homogeneous 
deformations in the absence of defects (such as surfaces). For the ID model Q this 
yields 

POO 

E ch (y) := / W{y') dx, for y G W£>\0, oo), (3) 
or equivalently, written in terms of the displacement gradient u = y' — 1, 

poo 

E ch (y)= W(l + u)dx, for u E L 1 (0, oo), 
Jo 

where W(r) = 0(r) + 0(2r) was already defined above. 

We consider a Pi finite element discretisation of the Cauchy-Born model. Let Xh '■= 
{X ,Xi, . . . } C N be a strictly increasing sequence of grid points with X Q = 0, and 
let hj := Xj + i — Xj. A Pi discretisation of y corresponds to a P discretisation of 
the displacement gradient u, hence we define for (Uj)'jL Q C R, where Uj denotes the 
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Figure 2. Visualisation of ([5]): the bond at the bottom of the graph 
is counted half in the Cauchy-Born model, even though it does not exist 
in the atomistic model, hence it gives a contribution — |0(2y'(O)) to the 
surface energy. 



displacement gradient in the element (X,-, Xj+i), 

oo 
3=0 

The Cauchy-Born approximation commits an error at the crystal surface, which 
the surface Cauchy-Born (SCB) approximation aims to rectify. The idea of the SCB 
method (in our ID setting) is to define 

POO 

E scb (y):= / W(y')dx + 1 (y'(0)), (4) 



and choose 7 such that the energy is exact under homogeneous deformations, which 
yields the formula 

7 (F) := -|0(2F); (5) 

see also Figure [2] Converting to the displacement gradient coordinate discretised by 
the Po finite element method we obtain 

E*>{V) :=E?(U) +1 (l + U ). 



Proposition 2. Ef^ and hence i£| are well-defined and twice Frechet differentiable 
in the weighted space i\(Xh) '■— {V = (Vj)'jL } equipped with the norm 

00 
3=0 

The first and second variations of E^ h are given by 

00 

(5E s h ch (U), V)=J2 W(l + UfiVi + 7 '(1 + Uo)V , 
3=0 

00 

(6 2 E s h cb (U) V, W}=J2 h 3W"(l + U 3 )V 3 W 3 + 7 "(1 + U )V W . 
3=0 
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2.3. Analysis of the linearized models. The parameter r for the Morse potential 
was chosen so that 1 is the minimizer of the Cauchy-Born stored energy function, which 
implies that 

t/f:=0, forj = 0,l,... (6) 

is the ground state of E^ 3 . More generally, u ch := (0)^ gives the bulk ground state of 
the crystal described by the model Q. We now consider linearisations of Eff b and E 3, 
about the Cauchy-Born state: 5E(0) + S 2 E(0)u = 0, where E G {E a ,E h h ch }. 
From Proposition |2] we obtain the linearised optimality condition for Eff h , 

7 / (l) + (W , (l)+7"(l))f/o = 0, and 

hjW"(l)Uj = for j = 1,2,..., 

which gives the linearised surface Cauchy-Born solution 

From Proposition [T] we obtain the linearised optimality condition for the atomistic 
model _E a , 

0'(1) + 0"(l)u o + (f>'(2) + (j>"{2){u Q + Ul ) = 0, 
0'(1) + (j)"{l)u 3 + 20'(2) + 0"(2)(u i _i + 2 Uj + u j+1 ) = 0, j > 1, 

which, using the fact that 0'(1) + 20'(2) = W (1) = can be rewritten in the form 

[0"(l) + 0"(2)K + 0"(2K = / (2), 
0"(2H-i + W'{1) + 20"(2)H + 0"(2)«, +1 = 0, £ > 1. 
This finite difference equation can be easily solved explicitly, which yields the solution 



m« := — — — - - }\ . — , where A = — tAA. fg\ 

1 0"(1) + 0"(2)(1 + A)' K) 

is the unique solution in (0, 1) of the characteristic equation 

0"(2)A 2 + [0"(1) + 20"(2)]A + <f>"{2) = 0. 

Since the expressions for ^ and ^ are somewhat bulky we expand them in the 
stiffness parameter a. The rationale for expanding in this parameter is that all models 
should coincide in the limit a — > 00. We hope, however, that our results will also yield 
useful predictions for moderate a. The elementary proof is postponed to the appendix. 

Proposition 3. Asymptotically as a — >■ we have the expansions 

Ut = S [1 - (1 + l)e" a + 0{e^)] , and (9) 
ul= ^[l-4e- a + 0(e- 2a )}. (10) 
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Remark 1. The asymptotic expansions (|9]) and (10) justify a posteriori the lineari- 
sation since they show that the displacements from the Cauchy-Born state are indeed 
small in the limit as a — > oo. □ 

2.4. Error estimates. We first note that each Po function U = (C/j)?l can be under- 
stood as a lattice function u = (i^)^L through the interpolation 

U£ = Uj for £ = Xj,. . ., X j+1 - 1, j e N. 

With this interpolation we obtain u ch = and w scb from the linearized CB and SCB 
solutions U ch and U scb , given in ([7]). 

We are interested in the improvement the SCB model gives over the pure Cauchy- 
Born model, that is, we wish to measure the relative errors 

||ti scb — M a L P ||li scb — U & \\iv 

Err„ := -. — r n — = n — n . 

||m cd — u a \\tP \\ ua \W 

Of particular interest are the uniform error Err^ and the error in the energy- norm Err2. 
We shall consider two separate cases: ho > 1 and ho = 1. 

Proposition 4 (Strain error). Let p G [l,oo] and h > 1, then 

Err p = C p + 0(e- a ), (11) 
where | < C p < 2. If ho = 1, then 

Err p = 2 1/p e~ a + 0(e- 2a ). (12) 



Proof. We consider the case h = 1 first. In that case (10) gives us 

1/p , oo \ 1/p 



E 



uf h - u a A p 



E 



it 



a \p 



A«S(1-A')' 



-i/p 



and similarly, ||w a ||£p = Uq(1 — \ p ) 1 ^ p . Using the asymptotic expansions (18) for A it is 
straightforward to show that 

(1 - \ p )- 1/p = 1 + 0(A) = 1 + C(e- a ); 



hence employing also (10) we obtain 



u 



a 



+ 0(^), and £ 



\uf b -uT 



i=i 



1 /p „-2a 



a 



+ 0{^). (13) 



For I = 0, since /iq = 1> we have 



k cb -<l 



r[l-3e- + 0(e- 



2a 



^[l-4e-« + 0(e- 2a )] = ^ + 0(^). 
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Combined with (13) this gives 
Err„ 



\u 



-.scb 1 1 
U \\f P 



\u a h P 



g \ g I 

g \ g I 



2 1 /P e - a + 0(e~ 2a ), 



which concludes the proof of (12). 

In the case h > 1 the convenient cancellation of first-order terms in u s Q cb — does 
not occur. Instead, using (13) we obtain 



\u a -u sch \ 



a 



Xi-l 
ho I / -i I ho I 

i=i 



i/p 



This immediately gives (11). 



□ 



We see from (11) that if we use a coarse finite element mesh up to the boundary, 
then the error in the displacement gradient will be typically of the order 50% or more. 
By contrast, if we refine the finite element mesh to atomistic precision at the boundary 
then the relative error is exponentially small in the stiffness parameter a. 

The quantity Err p measures the error in a pointwise sense. However, in some cases 
we are only interested in correctly reproducing certain macroscopic quantities such as 
the mean strain error 

Er= « b - % a ) 



Err : = 



2^=o u e 



Note that, up to higher order terms, this error also bounds the error in the displacements 
at the finite element nodes, which we observed in Figure [T] to be much smaller than the 
strain error. 

In the following result we confirm that, indeed, the mean strain error is an order of 
magnitude smaller than the pointwise strain error. 



Proposition 5 (Mean strain error). 

error satisfies 



Asymptotically as a — » oo, the mean strain 



Err = 2(1 - f )e~ a + 0(e 



-2g\ 



(14) 



Proof. We first compute the mean strains in the atomistic and the SCB models. For 
the atomistic model we have 



oo 

u" := > u 



■a . „,a . ^0 



^ 1 - A 

1=0 



Since (1 - A)" 1 = 1 + e~ a + 0{e~ 2a ) we obtain 



u & = ^ [(1 - 4e- Q )(l + e~ a ) + 0(e~ 2a )] = [l - 3e~ a + 0(e~ 2a )] . 
For the SCB model, we have 

oo 

-scb \ v u rrscb u rrscb e~ a |~i 1 1 i 2 \ „— a i p\( — 2a 



J2 hVf = h Ut = ^ [1 - (1 + £) c - + 0(e 

3=0 
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Figure 3. Relative error in the W /1 ' 2 -seminorm of the ID nonlinear SCB 
model for varying stiffness parameter a and two types of finite element 



grids; cf. Section 2.5 



and hence the error is given by 

T scb 



U 



U 



2(l-£>^ + C(^) 



This immediately implies (14) 



□ 



Remark 2. Since E 3, and E sch are Frechet differentiable in suitable function spaces 
it should be possible, using nonlinear analysis techniques such as the inverse function 
theorem, to extend the results from the linearized model problem to the fully nonlinear 
problem, provided that the stiffness parameter a is sufficiently large. Techniques of this 
kind have been used, for example, in [32]. □ 

2.5. Numerical results. We confirm through numerical experiments that the results 
of Propositions [4] and [5] are still valid in the nonlinear setting. In these experiments 
we choose tq — 1 instead of ([2]), choose a finite chain with 31 atoms, and let a vary 
between 2 and 7. For experiments with ho = 5 the gridpoints for the Cauchy-Born 
and SCB models are chosen as X = (0, 5, 10, ... , 30). For experiments with ho = 1, the 
gridpoints are chosen as X = (0, 1, 5, ... , 25, 29, 30). 

The results of the experiments are displayed in Figures [3] and |4j All results except 
for the relative error in the mean strain with ho — 1 confirm our analytical results in 
the linearized case. We have, at present, no explanation why the mean strain error Err 
with ho = 1 is of the order 0(e~ 3a ) instead of 0(e~ 2a ). A finer asymptotic analysis in 
the linearized case would in fact give the expansion Err = 2e~ 2a + 0(e~ 3a ). 
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Figure 4. Relative error in the mean strain of the ID nonlinear SCB 
model for varying stiffness parameter a and two types of finite element 



grids; cf. Section 2.5 



3. Numerical Results in 2D 

In this section we investigate numerically, to what extent the ID results might extend 
to the 2D setting. We will formulate a problem in a semi-infinite strip, where we expect 
relaxation only in the normal direction to the surfaces. Hence the problem reduces to 
a ID problem for the displacements in that direction. The ID analysis can be applied 
to this case with only minor changes, and we therefore expect the same behaviour as 
in the ID case. This is fully confirmed by the results of our numerical experiment. 

3.1. Formulation of the SCB method. In 2D one expects (this is rigorously proven 
only for large stiffness parameter a [31]) that the ground-state under Morse poten- 
tial interaction is the triangular lattice. Hence we choose as the atomistic reference 
configuration a subset A C AZ 2 , where 



1/2' 
x/3/2 



For future reference, we define a± := (1,0), 02 := (1/2, a/3/2) and 03 := (—1/2, -y/3/2), 
which are the directions of nearest-neighbour bonds. 
Specifically, we choose Ni,N 2 G N and define 

A : = {A(ri!,n 2 ) T G Z 2 | 1 < m < N u < n 2 < N 2 }, 

as the periodic cell of the semi-infinite strip A* := {A(n 1; n 2 ) T E Z 2 | < n 2 < N 2 }] cf. 
Figure [5j The corresponding continuous domain is Vt := A((0, A^] x (0, A^]). 
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Figure 5. Computational domain used in the numerical experiment 
described in Section [3} The small disks denote the set A; the dotted grid 
is the micro-triangulation T a ; the the large black disks denote the finite 
element nodes; the large white discs denote finite element nodes that are 
periodically repeated; the black lines denote the macro-triangulation T h . 



An admissible deformed configuration is a map y : A* — > IR 2 , which is periodic in 
the ai-direction, that is, y(£ + A^ai) = y(£) + Nidi. 

For simplicity we consider only second-neighbour interactions (measured in hopping 
distance). For each £ e A let jVg := {i] 6 A* | \r) — £| < 2} denote the interaction 
neighbourhood of £, then the potential energy of a deformed configuration is given by 

where <ft is again the Morse potential. 

To evaluate the deformation gradient dy of a discrete deformation y, we note that 
A* has a natural triangulation 71 (see Figure and identify y with its continuous 
piecewise affine interpolant in Pi(%,]M. 2 ). 

Let Th be a coarse triangulation of Q (which can be repeated periodically) and let 
Pi(7/ l ; M 2 ) denote the space of continuous and piecewise affine deformations of Q, such 
that yh(x+Niai) = yh,(x)+Niai, then the SCB energy of a deformation y^ e Pi (7^,; M 2 ) 
is given by 

E scb (y h )= [ W(dy h )dx+ [ 7 (9y h , z/) dar, 
Jn Jr 

where T C <9f2 denotes the free boundary, that is the portion of the boundary with 
normal v = ±(0,1), W is the Cauchy-Born stored energy function and 7 the SCB 
surface energy function, which are defined as follows: 

• If we denote by A/"cb the interaction neighbourhood of the origin in the infinite 
lattice AZ, 2 (see Figure [6^a)), then the Cauchy-Born stored energy function is 
given by 
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(a) 




Figure 6. (a) Third interaction neighbourhood, (b) Construction of 
7: Bonds A, B are underestimated by the Cauchy-Born approximation 
(counted only half), while the bonds C, D, E are overestimated (they do 
not exist in the atomistic model but are counted half in the Cauchy-Born 
model). 



• To define 7, we assume throughout that all surfaces of f2 are aligned with one 
of the three directions ai, a%, or a 3 , that is, u _L aj =: v L . Then the requirement 
that the SCB energy is exact under homogeneous deformations, in domains 
without corners, yields the expression 

7(F, 7 ) = ^(|F^|) + i0(2|F^|) 

- |0(V3|Fzy|) - §0(2|FQ 12 z/|) - §0(2|FQ?>|), 

where Q12 denotes a rotation through arclength 27r/12; see Figure [6](b) for an 
illustration. A rigorous proof of this formula follows immediately from Shapeev's 
bond density lemma [29] . 

3.2. Numerical results. In the numerical experiments we consider two types of finite 
element grids: a uniform grid with spacing h = ho = 5 (cf. Figure [Sta)), and a grid 
with an additional layer of elements at the free boundary, atomic spacing ho = 1 in 
the normal direction and uniform spacing h — 5 in the tangential direction (cf. Figure 
gb)). We will again measure the following relative errors: 

Err . /: _ \\dyl ch -dy*\\L> and I^hr dy»)dx 



fh - dy a \\L2 

where y a ,y s ^ h , and y^ 3 denote the minimizers of, respectively, E a ,E sch , and E sch with 
7 = 0. That is, Err 2 and Err measure the improvement of SCB over the pure Cauchy- 
Born model. 

The numerical results are displayed in Figures [7] and [8] Although the numerical 
results do not as clearly display the predicted convergence rates, they do seem to ap- 
proach these rates for increasing values of a. What is again clear is that the average 
strain has a much higher accuracy than the strain field, and that the additional mesh 
layer also substantially improves the accuracy of the method. We also note that we now 
observe essentially the predicted rate e~ 2a for the mean-strain error of the enhanced 
SCB model, instead of the unexpected rate e~ 3a . 



L( d y c h ~dy a ) dx 
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Figure 7. Relative error in the H^ 1,2 -seminorm of the 2D SCB model 
in the flat interface example described in Section |3j for varying stiffness 
parameter a and two types of finite element grids. 



Mean Strain Error 




Figure 8. Relative error for the mean strain of the 2D SCB model 
applied to the flat interface example described in Section [3j for varying 
stiffness parameter a and two types of finite element grids. 



Conclusion 



We presented an error analysis of the SCB method in the case where the dominant 
effect is surface relaxation in the normal direction. Our main results are: 1. We 
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showed that the "correct" approximation parameter is the stiffness of the interaction 
potential. 2. We showed that the mean strain (which is an important quantity of 
interest) has a much lower error than the strain field. 3. We showed that adding a 
single mesh layer at the free boundary with atomic spacing in the normal direction 
yields a substantial improvement to the accuracy of the SCB method with minimal 
increase in the computational cost. 

We also performed numerical experiments for domains with corners, which remain 
inconclusive so far. At corners there is an interplay between the normal stress and tan- 
gential stress of adjacent edges, which creates additional elastic fields. A finer analysis 
of this case is still required. In particular, it would be interesting to understand whether 
normal or tangential forces dominate the bahaviour of the system in that case. 

Appendix A. Proofs 

Proof of Propositions^ and\^ For each i G N we have 

0(1 + u t ) + 0(2 + u t + u e+1 ) = 0(1) + 0'(1H + §0"(0? } )M 2 

+ 0(2) + <f>'(2)(ut + u /+ i)^"(0f)k + ui+i\ 2 . 

where (9g — j) G i 1 by Taylor's theorem. Since 0(1) + 0(2) = 0, summing over £ G N 
and noting that the first-order terms cancel, yields 

E*(y) < C\\u\\%. 

Note that this seemingly requires only that uGf 2 , however, the series converges abso- 
lutely only if u G i 1 . 

Repeating the argument for a perturbation from a general state E a (u + v) shows the 
Frechet differentiability of E a . 

The same argument can be applied to prove Proposition [2j □ 

Proof of Proposition^ Inserting the definition of ro from ^ into 0"(1) yields 

0"(1) = We- 2a{1 - ro) - 2a 2 e~ a{1 - ro) 



l + 2e~ 2a J Vl + 2e- 2Q 



Expanding 



we obtain 



1 + 2e~ a 
1 + 2e~ 2a 



1 + 2e~ a + 0(e 



-2a\ 



0"(1) = 4a 2 (l + 4e" a ) - 2a 2 (l + 2e~ a ) + £>(aV 2Q ) 

= 2a 2 + 12a 2 e~ a + 0{a 2 e~ 2a ). (15) 
Similar calculations yield the expansions 

0'(2) = 2ae~ a + 2ae~ 2a + 0{ae~ 3a ), and (16) 
0"(2) = - 2a 2 e~ a + 0(a 2 e~ Za ). (17) 
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Writing out UQ Ch in terms of the Morse potential, and using the fact that 2 < 4 
2/ h < 4, which ensures that <f>"(l) + (4 - 2/h )4>"(2) > W"(l) > 0, we obtain 

0'(2) 0'(2) 1 



h U 



<j>»{\) + (4 - 2/M0"(2) 1 + (4 - 2/h )^ 

m 
0"(i) 



1 _ f4_l)ffl + offffl^ 



h J </>"(!) 



Inserting the expansions (15) to (17) gives (pi). 
To prove (10) we first expand A in terms of /3 :- 



and then in terms of e Q , 



A = 2^(^1 + 4/3-1-2/?) 

= i(l + 1(4/3) - |(4/3) 2 + ^(4/3) 3 + 0(/3 4 ) -1-2/?) 
= - /9 + 2/3 2 + C(/3 3 ) = e" a - 4e~ 2a + £>( e - 3a ). 
Inserting this result into pi) and a brief computation yield 



0'/ + 0' 2 '(1 + A) 



\ a V a I 



Since = u%\ the result (1101) follows easily. 



□ 
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